Ray dynamics in a long-range acoustic propagation experiment 
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A ray-based wavefield description is employed in the interpretation of measurements made during 
the November 1994 Acoustic Engineering Test (the AET experiment). In this experiment phase- 
coded pulse-like signals with 75 Hz center frequency and 37.5 Hz bandwidth were transmitted near 
the sound channel axis in the eastern North Pacific Ocean. The resulting acoustic signals were 
recorded on a moored vertical receiving array at a range of 3252 km. In our analysis both mesoscale 
and internal-wave-induced sound speed perturbations are taken into account. Much of this analysis 
exploits results that relate to the subject of ray chaos; these results follow from the Hamiltonian 
structure of the ray equations. We present evidence that all of the important features of the measured 
AET wavefields, including their stability, are consistent with a ray-based wavefield description in 
which ray trajectories are predominantly chaotic. 
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I. INTRODUCTION 

The chaotic dynamics of ray trajectories in ocean 
acoustics have been explored in a number of recent pub- 
lications (Q-H2), including two recent reviews [2lL l22fl. 
In these publications a variety of theoretical results are 
presented and illustrated, mostly using idealized mod- 
els of the ocean sound channel. In contrast, in the 
present study a ray-based wavefield description is em- 
ployed, in conjunction with a realistic environmental 
model, to interpret a set of underwater acoustic mea- 
surements. For this purpose we use both oceanographic 
and acoustic measurements collected during the Novem- 
ber 1994 Acoustic Engineering Test (the AET experi- 
ment). In this experiment phase-coded pulse-like signals 
with 75 Hz center frequency and 37.5 Hz bandwidth were 
transmitted near the sound channel axis in the eastern 
North Pacific Ocean. The resulting acoustic signals were 
recorded on a moored vertical receiving array at a range 
of 3252 km. We present evidence that all of the impor- 
tant features of the measured AET wavefields are consis- 
tent with a ray-based wavefield description in which ray 
trajectories are predominantly chaotic. 

The AET measurements have previously been ana- 
lyzed in Refs. [2j|-[2f|. In the analysis of Worcester 
et al. [2^| it was shown that the early AET arrivals 
could be temporally resolved, were stable over the du- 
ration of the experiment and could be identified with 
specific ray paths. These features of the AET observa- 
tions are consistent with those of other long-range data 
sets [Mill- Tn e AET arrivals in the last 1 to 2 sec (the 
arrival finale) could not be temporally resolved or iden- 
tified with specific ray paths, however. Also, significant 
vertical scattering of acoustic energy was observed in the 
arrival finale. These features of the AET arrival finale 



are consistent with measurements and simulations (both 
ray- and parabolic-equation-based) at 250 Hz and 1000 
km range [2^ | -|29] | and at 133 Hz and 3700 km range [30| . 

The analysis of Colosi et al. [2J| focused on statistical 
properties of the early resolved AET arrivals. Measure- 
ments of time spread, travel time variance and proba- 
bility density functions (PDFs) of peak intensity were 
presented and compared to theoretical predictions based 
on a path integral formulation as described in Ref. |3l| . 
Travel time variance was well predicted by the theory, 
but time spreads were two to three orders of magnitude 
smaller than theoretical predictions, and peak intensity 
PDFs were close to lognormal, in marked contrast to the 
predicted exponential PDF. The surprising conclusion of 
this analysis was that the measured AET pulse statis- 
tics suggested propagation in or near the unsaturated 
regime (characterized physically by the absence of mi- 
cromultipaths and mathematically by use of a perturba- 
tion analysis based on the Rytov approximation), while 
the theory predicted propagation in the saturated regime 
(characterized by the presence of a large number of mi- 
cromultipaths whose phases are random). 

In Ref. (2j| it was shown that in the finale region of 
the measured AET wavefields, where no timefronts are 
resolved, the intensity PDF is close to the fully saturated 
exponential distribution. This result is not unexpected 
because the finale region is characterized by strong scat- 
tering of both rays [l4| an( A modes j^J ■ 

The work reported here seeks to elucidate both the 
physics underlying the AET measurements, and the 
causes of the successes and failures of the analyses re- 
ported in Refs. [2j|-[2{|. We show that a ray-based anal- 
ysis, in which ideas associated with ray chaos play an 
important role, can account for the stability of the early 
arrivals, their small time spreads, the associated near- 
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lognormal PDF for intensity peaks, the vertical scatter- 
ing of acoustic energy in the reception finale, and the 
near-exponential intensity PDF in the reception finale. 
Part of the success of this interpretation comes about 
due to the presence of a mixture of stable and unstable 
ray trajectories. Also, an explanation is given for differ- 
ing intensity statistics in the early and late portions of 
the measured wavefields, and the fairly rapid transition 
between these regimes. 

The remainder of this paper is organized as follows. 
In the next section we describe the AET environment 
and the most important features of the measured and 
simulated wavefields. In the simulated wavefields both 
measured mesoscale and simulated internal- wave-induced 
sound speed perturbations are taken into account. It is 
shown that even in the absence of internal waves, the late 
arrivals are not resolved and the associated ray paths are 
chaotic. In section II we present results that relate to 
the structure of the early portion of the measured time- 
front and its stability. The micromultipathing process is 
discussed, and a quantitative explanation for the remark- 
ably small time spreads of the early ray arrivals is given. 
Section III is concerned with intensity statistics. An ex- 
planation is given for the cause of the different wavefield 
intensity statistics in the early and late portions of the 
arrival pattern. In the final section we summarize and 
discuss our results. 



assume propagation in the vertical plane with cartesian 
coordinates z and r representing depth and range, re- 
spectively. Ray methods can be introduced when the 
wavelength 2irc/a of all waves of interest is small com- 
pared to the length scales that characterize variations in 
c. Substituting the geometric ansatz, 



(2) 



representing a sum of locally plane waves, into the 
Helmholtz equation gives, after collecting terms in de- 
scending powers of a, the eikonal equation, 



(VT) 2 = <T 2 , 



and the transport equation, 

V(A 2 VT) = 0. 



(3) 



(4) 



For notational simplicity, the subscripts on A and T have 
been dropped in © and J3J. 

The solution to the eikonal equation can be reduced 
to solving a set of ray equations. If a one-way formulation 
(see, e.g., Ref. |22j) is invoked, with r playing the role of 
the independent variable, the ray equations are 



dz 
dr 



dH dp 
dp ' dr 



dH 

dz ' 



(5) 



II. MEASURED AND SIMULATED 
WAVEFIELDS 

Figure 1 shows a measured AET wavefield in the time- 
depth plane and ray-based simulations of such a wave- 
field with and without internal waves. Figure 2 is a com- 
plementary display that shows measured and simulated 
wavefields after plane wave beamforming. In this section 
we describe the most important qualitative features of 
measured and simulated AET wavefields. We begin by 
reviewing fundamental ray theoretical results so that our 
ray simulations can be fully understood. Also, this mate- 
rial provides a foundation for some of the later discussion. 
We then briefly describe the AET environment and the 
signal processing steps that were performed to produce 
the measured wavefield shown in Fig. 1. Next, our treat- 
ment of internal- wave- induced sound speed perturbations 
is briefly described. Finally, we return to describing and 
comparing qualitative features of the measured and sim- 
ulated wavefields. 

Fixed-frequency (cw) acoustic wavefields satisfy the 
Helmholtz equation, 



V 2 u + cr 2 c _2 (z,r)u = 0, 



(1) 



where a — 2irf is the angular frequency of the wave- 
field and c(z, r) is the sound speed. (The extension to 
transient wavefields is straightforward using Fourier syn- 
thesis. It is not necessary to consider this complication 
to present the most important results needed below.) We 



and 



where 



— - L - —-H 

dr ^ dr ' 



H(p,z,r) = - \jc 2 (z,r) 



V 



(6) 



(7) 



These equations constitute a nonautonomous Hamilto- 
nian system with one degree of freedom; z and p are 
canonically conjugate position and momentum variables, 
r is the time-like variable, H is the Hamiltonian and the 
travel time T corresponds to the classical action. Note 
also that p = p z — dT/dz and p r = dT/dr = —H are 
the z- and r-components, respectively, of the ray slow- 
ness vector. The ray angle relative to the horizontal ip 
satisfies dz/dr = tamp, or, equivalently, cp = simp. For a 
point source rays can be labeled by their initial slowness 
Pq. The solution to the ray equations is then z(r;po), 
p(r;po), T(r;po). The terms in (0) correspond to eigen- 
rays; these satisfy z(rn;po) = zr where the depth and 
range of the receiver are zr and rfj, respectively. The 
transport equation can be reduced to a statement of 
constancy of energy flux in ray tubes. The solution, for 
the J-th eigenray, can be written 



A j (z,r)=A \H/q 21 \) /2 



(8) 



where H , the matrix element q 2 i and the Maslov index 
/x are evaluated at the endpoint of the ray, and Aq is a 
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constant. The stability matrix 



921 922 



quantifies ray spreading, 
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Elements of this matrix evolve according to 
d 



dr 



Q = KQ 



where Q at r = is the identity matrix, and 



K = 



d 2 H d'H 

dzdp dz 2 

d 2 H d 2 H 

dp' 2 dzdp 



(9) 



(10) 



(11) 



(12) 



At caustics q 2 i vanishes and, for waves propagating in 
two space dimensions, the Maslov index advances by one 
unit. 

In ocean environments with realistic range- 
dependence, including the AET environment, ray 
trajectories are predominantly chaotic. Chaotic rays 
diverge from neighboring rays at an exponential rate, on 
average, and are characterized by a positive Lyapunov 
exponent, 



lim 



1 lim 
r 2?(0)- 



>0 T>(0)J 



(13) 



Here T>(r) is a measure of the separation between rays at 
range r; suitable choices for T> are any of the four matrix 
elements of Q or the trace of Q. The chaotic motion of ray 
trajectories leads to some limitations on predictability. 
This will be discussed in more detail below. Additional 
background material relating to this topic can be found 
in 0, and references therein. 

We digress now to briefly describe the AET experi- 
ment and relevant signal processing. In this experiment 
a source, suspended at 652 m depth from the floating 
instrument platform R/P FLIP in deep water west of 
San Diego, transmitted sequences of a phase-coded signal 
whose autocorrelation is pulse-like (with a pulse width 
of approximately 27 ms). The source center frequency 
was 75 Hz and the bandwidth was approximately 37.5 
Hz. The resulting acoustic signals were recorded east of 
Hawaii at a range of 3252.38 km on a moored 20 element 
receiving array between the depths of 900 m and 1600 m. 
After correcting for the motion of the source and receiv- 
ing array, including removal of associated Doppler shifts, 
the received acoustic signals were cross-correlated with a 
replica of the transmitted signal (to achieve the desired 
pulse compression) and complex demodulated. To im- 
prove the signal-to-noise ratio, 28 consecutive processed 
receptions, extending over a total duration of 12.7 min- 
utes, were coherently added. The duration over which 



this coherent processing yields improved signal-to-noise 
is limited by ocean fluctuations, mostly due to internal 
waves. The statistics of the acoustic fluctuations were 
shown in Ref. j2f| not to be adversely affected by the 12.7 
minute coherent integration. Nearly concurrent with the 
week-long experiment, temperature profiles in the upper 
700 m were measured with XBT's every 30 km along 
the transmission path. An objective mapping technique 
was applied [2^ to these measurements to produce a 
map of the background sound speed structure (including 
mesoscale structure) along the transmission path. The 
sampling interval in range of the resulting sound speed 
field is 30 km so that no structure with horizontal wave- 
lengths less than 60 km is present in this field. More 
details on the AET experiment, environment and pro- 
cessing of the acoustic signals can be found in Ref. [23| . 

Internal-wave-induced sound speed perturbations are 
taken into account in most of our simulations. These 
are assumed to satisfy the relationship 5c = {dc/dz) p ( > 
and the statistics of £, the internal-wave-induced ver- 
tical displacement of a water parcel, were assumed to 
be described by the empirical Garrett-Munk spectrum 
(see, e.g., Ref. 1^3 )• The potential sound speed gradi- 
ent (dc/dz) p , the buoyancy frequency N, and the sound 
speed c were estimated directly from hydrographic mea- 
surements using a procedure similar to that described in 
the appendix of Ref. [3(| ; it was found that a good ap- 
proximation for the AET environment is (dc/dz) p /c « 
(1.25 s 2 /m)iV 2 . The vertical displacement C,(r,z) was 
computed using equation 19 of Ref. [34j with the vari- 
able x replaced by r, and y = t = 0. Physically this 
corresponds to a frozen vertical slice of the internal wave 
field that includes the influence of transversely propa- 
gating internal wave modes. A Fourier method, with 
Ak x = 2-7r/3276.8 km and < k. x < 2ir/l km, was used 
to numerically generate sound speed perturbation fields. 
(Ray calculations were also performed using an internal 
wave field in which the hard upper bound on k x was 
replaced by an exponential damping of the large k x en- 
ergy. The qualitative ray behavior described below did 
not depend on the manner in which the large k x cut-off 
was treated.) A mode number cut-off of j ma x = 30 was 
used in all the simulations shown. Contributions from ap- 
proximately 2 30 (2 15 / 5 x 2 15 / 5 x 30 where the factors of 
1/5 are present because the assumed wavenumber cut-off 
was one-fifth of the assumed Nyquist wavenumber) inter- 
nal wave modes were included in our simulated internal- 
wave-induced sound speed perturbation fields. Measure- 
ments of sound speed variance made during the experi- 
ment (see Fig. 3) indicate that the internal wave energy 
strength parameter E was close to Egm, the nominal 
Garrett-Munk value. Simulations were performed using 
both E = 0.5 Egm and E = 1.0 E GM - All of the simu- 
lations shown in this paper use E = 0.5 Egm- The dif- 
ference between these simulations and those performed 
using E = 1.0 Egm will be discussed where appropriate. 

With the foregoing discussion as background, we are 
prepared to discuss our ray simulations in the AET en- 
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vironment. These are shown in Figs. 1, 2, 4 through 
11, and 14. The points plotted in Figs. 1, 2 and 4 were 
computed by integrating the coupled ray equations (0 
from the source to the range of the receiving array. Sev- 
eral hxed and variable step-size integration algorithms 
were tested. In the presence of internal waves, conver- 
gence using a hxed step fourth order Runge Kutta in- 
tegrator required that the step size not exceed approxi- 
mately one meter. All ray tracing calculations presented 
in this paper were performed using double precision arith- 
metic (64 bit floating point wordsize). The waveheld in- 
tensity is approximately proportional to the density of 
rays (dots) that are plotted in Figs. 1 and 2. A some- 
what more difficult calculation in which the contributions 
to the waveheld from many rays are coherently added will 
be discussed below. 

In the presence of internal waves, rays with launch an- 
gles between approximately ±5° form the diffuse hnale 
region of the arrival pattern arriving after approximately 
2196 sec. Steeper rays contribute to the earlier, mostly 
resolved arrivals. It is seen in Fig. 1 that our ray sim- 
ulation in the presence of internal waves underestimates 
the vertical spread of the hnale region. Simulations with 
a stronger internal wave held, E = 1.0 Egm, give a bet- 
ter fit to the vertical spread of energy (but a poorer fit 
to other wavefield features, as described below). Note, 
however, that diffractive effects, missing in our ray simu- 
lations, will contribute to enhanced broadening in depth 
of the arrival pattern. A striking feature of Figs. 1, 2 
and 4 is the contrast between the seemingly structure- 
less distribution of steep (large \p\) rays in (p, z) seen in 
Fig. 4 when internal waves are present and the relatively 
structured distribution of the same rays in the time-depth 
(Fig. 1) and time-angle (Fig. 2) plots. 

The cause of the complexity seen in Fig. 4 is ray 
chaos [^1 |22; most ray trajectories diverge from ini- 
tially infinitesimally perturbed rays at an exponential 
rate. 48,000 rays with launch angles between ±12° were 
traced to produce Figs. 1, 2 and 4 in the presence of in- 
ternal waves. This fan of rays is far too sparse to resolve 
what should be an unbroken smooth curve - a Lagrangian 
manifold - which does not intersect itself in phase space 
(Fig. 4). Under chaotic conditions the separation be- 
tween neighboring rays grows exponentially, on average, 
in range. The complexity of the Lagrangian manifold 
grows at the same exponential rate. The Lyapunov expo- 
nent is the recriprocal of the e- folding distance (see, e.g., 
|22|). Finite range numerical estimates of Lyapunov ex- 
ponents (hereafter referred to as stability exponents) are 
shown as a function of launch angle in Fig. 5. It is seen 
that in this environment the flat rays (\<fo\ ~ 5°) have sta- 
bility exponents of approximately 1/(100 km), while the 
steeper rays (6° $ \y>o\ $ 11°) have stability exponents of 
approximately 1/(300 km). It follows, for example, that 
the complexity of the flat ray portion of the Lagrangian 
manifold grows approximately like exp(r/100 km). 

Not surprisingly, stability exponents show some sen- 
sitivity to parameters that describe the internal wave 



field. Our simulations show that the stability exponents 
of rays with axial angles of approximately 5° or less ap- 
proximately double when j ma x increases from 30 to 100, 
while steeper rays show very little sensitivity to j max - 
In contrast, stability exponents of rays with axial angles 
from approximately 8° to 15° approximately double when 
kmax increases from 27r/1000 m to 27r/250 m, while the 
flatter rays show little sensitivity to k max . With these 
comments in mind, one should not attach too much sig- 
nificance to the details of Fig. 5. We note, however, that 
rays with axial angles of approximately 5° or less consis- 
tently have higher stability exponents than rays in the 5° 
to 10° band. 

Note that what appears to be temporally resolved ar- 
rivals in the measured wavefield correspond in our simu- 
lations to contributions from an exponentially large num- 
ber of ray paths. (Only a small fraction of the total num- 
ber are included in our simulations, however, because of 
the relative sparseness of our initial set of rays.) The ob- 
servation that the travel times of chaotic ray paths may 
cluster and be relatively stable was first made in Ref. [g] • 

Simmen et al. |14| have previously produced plots sim- 
ilar to our Figs. 1 and 4 for ray motion in a deep ocean 
model consisting of a background sound speed struc- 
ture very similar to ours on which internal- wave-induced 
sound speed perturbations were superimposed. Although 
our results are similar in many respects, it is noteworthy 
that the right panel of our Fig. 4 shows more chaotic be- 
havior than is present in the corresponding plot in Ref. 
[T^. This difference persists if our 3252.38 km range sim- 
ulations are replaced by simulations at the same range 
(1000 km) that was used in Ref. 0]. O ur ra y s > espe- 
cially the steeper rays, are more chaotic than those in Ref. 
[l4| . We believe that this difference is primarily due to 
the increased complexity of our internal wave field over 
that of Simmen et al. |l4j who included contributions 
from ten internal wave modes in their simulated internal 
wave fields. 

Figures 6 and 7 show plots of ray depth at the AET 
range vs. launch angle. In Fig. 6 the same rays (exclud- 
ing those with negative launch angles) that were used to 
produce Figs. 1, 2 and 4 are plotted. In Fig. 7 two small 
angular bands of the upper panel of Fig. 6 are blown up; 
the ray density is four times greater than that used in Fig. 
6. Fig. 6 strongly suggests that, even in the absence of 
internal waves, the near-axial rays are chaotic. Not sur- 
prisingly, in the presence of internal waves, steeper rays 
are also predominantly chaotic. Note, however, that Figs. 
5, 6 and 7 show that the near-axial rays are evidently 
more chaotic than the steeper rays. An explanation for 
this difference in behavior will be given in section III. 

In this section we have briefly described and compared 
the gross features of measured and simulated AET wave- 
fields. Ray trajectories in our simulated wavefields in the 
presence of internal waves are predominantly chaotic. In 
spite of this, Figs. 1 and 2 show that many features of 
our simulated wavefields are stable and in good quali- 
tative agreement with the observations. Predicted and 
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measured spreads of acoustic energy in time, depth and 
angle are generally in good agreement, both in the early 
and late portions of the arrival pattern. We shall not 
discuss further the spreads of energy in depth and angle, 
except to note that Figs. 1 and 2 show that these spreads 
can be accounted for using ray methods in the presence 
of realistic (including internal-wave-induced sound speed 
perturbations) ocean structure. In the two sections that 
follow we shall consider in more detail the time spreads 
of the early ray arrivals, and intensity statistics in both 
the early and late portions of the arrival pattern. 



III. CHAOS, MICROMULTIPATHS AND 
TIMEFRONT STABILITY 

In this section we consider in some detail eigenrays, 
timefront stability and time spreads. We focus our at- 
tention on the early branches of the timefront, where 
measured time spreads were, surprisingly, only approxi- 
mately 2 ms, in sharp contrast to the theoretical predic- 
tion of approximately 1 sec [24j. We focus our attention 
on the early arrival time spreads, primarily because quan- 
titative estimates of the unresolved late arrivals are not 
available. 

Eigenrays at a fixed depth z correspond to the inter- 
sections of a horizontal line at depth z with the curve 
z(ip ); these intersections are the roots of the equation 
z — z(tp ) = 0. Although the discrete samples of z(ip ) 
plotted in Figs. 6 and 7 are too sparse to reveal what 
should be a smooth curve, it is evident the number of 
eigenrays at almost all depths is very large; because rays 
are predominantly chaotic, this number grows exponen- 
tially in range. In principle, eigenrays can be found us- 
ing a combination of interpolation and iteration, starting 
with a set of discrete samples of z(ip D ). In practice, this 
procedure reliably finds only those eigenrays in regions 
where z(ip Q ) has relatively little structure. These rays 
have the highest intensity and are the least chaotic. It is 
seen in Fig. 1 that the early (steeper) ray arrivals come 
in clusters with small travel time spreads. The set of 
eigenrays within such a cluster is often referred to as a 
set of micromultipaths. An incomplete set of micromulti- 
paths, found using the procedure just described, is shown 
in Fig. 8. An interesting and important feature of micro- 
multipaths in our simulations is that all of the rays that 
make up a set of micromultipaths have the same iden- 
tifier ±M. Here ± is the sign of the launch angle and 
M is the number of ray turning points (where p changes 
sign) following a ray. In Fig. 9 rays with two fixed val- 
ues of ray identifier (+137 and +151) are plotted in the 
time-depth plane. In both cases the points plotted are 
a subset of those plotted in Fig. 1. It is seen that the 
relatively steep +137 rays have a small time spread and 
form one of the branches of the timefront seen in Fig. 1 
(similar behavior was noted in Ref. 0]), while the rela- 
tively flat +151 rays have a much larger time spread and 
fall within the diffuse finale region of the arrival pattern 



seen in Fig. 1. 

An important (and surprising) feature of sets of micro- 
multipaths is that they are highly nonlocal in the sense 
that interspersed (i.e., having intermediate launch an- 
gles) among a group of micromultipaths are rays with M 
values that differ by several units. This behavior is illus- 
trated in Fig. 10. In view of the observation that, in the 
presence of internal waves in the AET environment, ray 
motion is strongly chaotic and the function M((p ) has 
local oscillations of several units, it is remarkable that 
the early portion of the timefront (see Fig. 1) is not de- 
stroyed. The nonlocality of a set of micromultipaths in 
the range-depth plane is seen in Fig. 8. Note that there 
are significant differences in the upper and lower turn- 
ing depths of the plotted micromultipaths and that these 
rays are spread in range by a significant fraction of a ray 
cycle distance. 

The result of performing a ray-based synthesis of what 
appears in the measurements to be an isolated arrival 
is shown in Fig. 11. This involves finding a complete 
set of micromultipaths, including their travel times and 
Maslov indices, and coherently summing their contribu- 
tions. Unfortunately, owing to the predominantly chaotic 
motion of ray trajectories in the AET environment in 
the presence of internal waves, it is extremely difficult 
to find a complete set of micromultipaths. Indeed, the 
constituent ray arrivals used in the Fig. 11 synthesis do 
not constitute a complete set of micromultipaths. This is 
less important than might be expected because standard 
eigenray finding techniques easily locate the strongest mi- 
cromultipaths; the missing micromultipaths in the Fig. 
11 synthesis are highly chaotic rays that have very small 
amplitudes. An interesting result of this synthesis is that 
the micromultipath-induced time spread in the synthe- 
sized pulse is only about 1 ms, which is in very good 
agreement with the AET measurements [3 HH- Note 
that the total time spread among the micromultipaths 
found - about 11 ms - is much greater than the spread 
of the synthesized pulse. This difference arises because 
the time spread of the dominant micromultipaths is much 
smaller than the total micromultipath time spread. 

Simulations (not shown) using E — 1.0 Eqm result 
in synthesized steep ray pulses that are spread in time 
much more than is shown in Fig. 11. With the stronger 
internal wave field, the total micromultipath time spread 
increases by approximately a factor of two - to about 
25 ms. More importantly, however, the dominant mi- 
cromultipaths have time spreads comparable to the to- 
tal micromultipath time spread, leading to synthesized 
pulses that are spread by 5 to 10 ms, rather than 1 to 
2 ms. Thus, measured early arrival time spreads are in 
better agreement with E = 0.5 Eqm simulations than 
with E = 1.0 Eqm simulations. The question of which 
value of the internal wave strength E in our simulations 
gives the best fit to the observations will be revisited in 
the next section. 

An interesting and unexpected feature of our simula- 
tions is that the Maslov index was consistently lower than 
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the number of turns M made by the same ray. The dif- 
ference was typically three units in the E = 0.5 Eqm 
simulations and four units in the E = 1.0 Eqm simula- 
tions. We have no explanation for this behavior. 

So far in this section we have seen that associated with 
the chaotic motion of ray trajectories is extensive micro- 
multipathing, and that the micromultipathing process is 
highly nonlocal. The many micromultipaths that add 
at the receiver to produce what appears to be a single 
arrival may sample the ocean very differently. Surpris- 
ingly, on the early arrival branches the highly nonlocal 
micromultipathing process causes only very small time 
spreads and does not lead to a mixing of ray identifiers. 
To complete this picture it is necessary to explain why 
the time spreads of the early (steep) ray arrivals are so 
small. We shall now provide such an explanation. Re- 
call first, however, that the measured time spread of the 
early arrivals is only approximately 2 ms J23, in marked 
contrast to the theoretical prediction 0, l3l| of approxi- 
mately 1 sec. Thus our quantitative explanation for the 
small early arrival time spreads serves both to give insight 
into the underlying physics, and provides an explanation 
for a critically importnat feature of the measured AET 
wavefields. 

Computing time spreads is conceptually straightfor- 
ward using ray methods. In a known environment one 
finds a large ensemble of rays with the same identifier 
that solve the same two point boundary value problem; 
the spread in travel times over many ensembles of such 
rays, each ensemble computed using in a different real- 
ization of the environment, is the desired quantity. The 
constraint that all of the computed rays have the same 
fixed endpoints complicates this calculation. In the fol- 
lowing we exploit an approximate form of the eigenray 
constraint that simplifies the calculation; the approxi- 
mate form of the eigenray constraint is first order accu- 
rate in a sense to be described below. 

Before describing the manner in which the eigenray 
constraint is imposed, it is useful to describe the physi- 
cal setting in which the constraint will be applied. Figure 
12 shows \H\ vs r following a moderately steep ray (axial 
angle approximately 11°) in a canonical environment |35l | 
on which the previously described (based on the mea- 
sured AET N(z) profile) internal-wave-induced sound 
speed perturbation field was superimposed. It is seen 
that H (r) consists of a sequence of essentially constant H 
segments, separated by near-step-like jumps. The jumps 
occur at the ray's upper turning points. A simple model 
that captures this behavior - the so-called apex approx- 
imation |3l| - assumes that the transition regions can 
be approximated as step functions. (The validity of the 
apex approximation is linked to the anisotropy and inho- 
mogeneity of oceanic internal waves; details of the back- 
ground sound speed profile are not important. For this 
reason we have chosen to illustrate this effect using a 
background canonical profile rather than the AET envi- 
ronment.) Figure 13 shows a schematic diagram of two 
rays. One ray has undergone two internal- wave-induced 



apex scattering events; the other ray is an unscattered 
ray with the same launch angle as the scattered ray. 

Our strategy to building the eigenray constraint into 
an estimate of travel time spreads can now be stated. We 
consider a scattered and an unperturbed (that is, in the 
absence of internal-wave-induced scattering events) ray 
with the same ray identifier; the rays shown schematically 
in Fig. 13 have identifier +3. In addition to having the 
same ray identifier, the scattered and unperturbed rays 
are assumed to start from the same point and end at 
the same depth, but they will generally end at different 
ranges. We assume that the scattered ray has travel time 
T r and is one of many scattered eigenrays with the same 
ray identifier at range r. We estimate the travel time T(r) 
of an unperturbed eigenray at range r whose ray identifier 
is the same as that of the scattered eigenray. It will be 
shown that 6T(r) = T r — T(r) vanishes to first order if 
the apex approximation is exploited. Because the same 
result applies to all of the scattered eigenrays, it follows 
that the time spread vanishes to first order, independent 
of r, if the apex approximation is exploited. Note that 
this result is not expected to apply, even approximately, 
to the late near-axial rays where the apex approximation 
is known to fail. 

The travel time T(r) of the eigenray in the unper- 
turbed ocean can, of course, be computed numerically 
but this provides little help in deriving an analytical es- 
timate of ST. Instead, consider a ray in the unperturbed 
ocean with the same launch angle as one of the scattered 
eigenrays, as shown in Fig. 13. In general, the range 
of this ray at the receiver depth, after making the ap- 
propriate number of turns, will be r Q ^ r. But T(r) 
can be estimated from T(r ). This follows from the ob- 
servation that ray travel time (assuming a point source 
initial condition, for instance) is a continuously differen- 
tible function of z and r with VT = p where p is the ray 
slowness vector. Because of this property T(z, r) can be 
expanded in a Taylor series. If we consider rays with the 
same ray identifier and fix the final ray depth to coincide 
with the receiver depth, then 

T{r)^T(r )+ Pr {r ){r-r ) (14) 

where p r = — H is the r-component of the ray slow- 
ness vector. (More generally, T(z,r) consists of multiple 
smooth sheets that are joined at cusped ridges; in spite 
of this complication, the property of continuous differen- 
tiability is maintained, even in a highly structured ocean. 
For our purposes, we need only consider one sheet of this 
multi-valued function. Note also that our application of 
(1141) in the background ocean makes this expression par- 
ticularly simple to use.) It follows from l|14H that 

5T(r) =T r - T{r) « T r - T(r ) - Pr (r )(r - r Q ) (15) 

is, correct to first order, the travel time difference at 
range r between eigenrays with and without internal 
waves. 

To evaluate 5T(r) we shall assume that the backgound 
environment is range-independent. Also, consistent with 
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the use of the apex approximation, we may treat the 
perturbed environment as piecewise range-independent. 
Within each range-independent segment it is advanta- 
geous to make use of the action-angle variable formalism 
(see. e.g., [2^ or |3(j). In terms of the action-angle vari- 
ables (I, 0), H = H(I) and the ray equations are 



dO dH 
d-r = W= W(7) ' 



dl 
dr 



dH 
~d0 



= 0, 



and 



(16) 



(17) 



(18) 



where the physical interpretation of H as — p r is main- 
tained. The angle variable can be defined to be zero 
at the upper turning depth of a ray and increases by 
2ir over a ray cycle (double loop). It follows that the 
frequency oj(I) = 2tt/R(I) where R(I) is a ray cycle 
distance. Integrating the ray equations over a com- 
plete ray cycle then gives R(I) = 2tt/u(I) and T(I) = 
2w(I — H(I)/oj(I)) with / constant following each ray. 
In the apex approximation / jumps discontinuously at 
each ray's upper turning depth; for the perturbed ray 
R(I + AI) k 2tt/uj(I) - 2ttAIuj / {I)/(lu{I)) 2 and T(I + 
AI) « 2w(I - H{I)/uj{I)) + 2ttAIH(I)oj'(I)/(uj(I)) 2 
where U)'(I) — du>(l)/dl. Note that, like these expres- 
sions, Eqs. 1141 andlTSlare first order accurate in AI. 

Consider again Eq. El an d Fig- 13 with I equal to 
lo, h and h in the left center and right ray segments, 
respectively. The left segment gives no contribution to 
ST as the ray has not yet been perturbed. In the center 
ray segment 



T r - T{r ) a 2tt(/i - J ) 



H(I )w'(I ) 
M/o)) 2 



and 



-Pr(r ){r~r ) sw H(I ) 



-2Tr(h - I a ) 



w'(Jo) 



(19) 



(20) 



These terms are seen to cancel. Note that if additional 
complete cycle ray segments are added to the center sec- 
tion of the ray, Eqs. I19landl2*0lare unchanged except that 
I\ — Iq = All is replaced by J27 =1 AT — I n — Iq; again 
the two terms cancel. In the final (incomplete cycle) ray 
segment the difference between the terms T r — T(r ) and 
Pr(r )(r — r ) can be shown to be 0((AI) 2 ); the terms do 
not exactly cancel because the ^-values of the perturbed 
and unperturbed rays are generally not identical at the 
receiver depth. Thus, to first order in AI, ST — 0, in- 
dependent of range, if the apex approximation is valid. 
This simple calculation provides an explanation of why, 
in spite of extensive ray chaos, the time spreads of the 
early AET ray arrivals are quite small. 



Several comments concerning the preceeding calcula- 
tion are noteworthy. First, we note that although it was 
assumed that the background sound speed structure is 
range-independent, the preceeding argument also holds 
in the presence of slow background range-dependence, 
i.e., with structure whose horizontal scales are large rela- 
tive to a typical ray double loop length. Adiabatic invari- 
ance in such environments guarantees that while H is not 
constant following rays between apex scattering events, 
/ is nearly constant. Second, after n random kicks I± — Iq 
in ijT§|) and 1(20) 1 is replaced by I n ~Io ~ y/n(AI) rms , and 
the magnitude of 1 (19) 1 under AET conditions is 0(1 sec). 
This is an example of a travel time spread estimate that 
fails to enforce the eigenray constraint. This calculation 
shows that the difference between travel time estimates 
that do and do not enforce the eigenray constraint can 
be quite significant. (The constrained estimate is refined 
below.) Third, it should be emphasized that Eqs. I|19|) 
and (|20|) apply (approximately) only to the steep rays be- 
cause the apex approximation applies only to the steep 
rays. We have not addressed the time spreads of near- 
axial rays. Note, however, that Fig. 9 shows that time 
spreads are larger for flatter rays. Fourth, the above cal- 
culation shows that, to lowest order in AI, there is no 
internal-wave-scattering-induced travel time bias if the 
apex approximation is strictly applied. And fifth, the ar- 
guments leading to Eqs. (|19fl and (|20|l apply whether the 
scattered rays are chaotic or not. 

A relaxed form of the apex approximation in which 
the action jump transition region has width A0 gives a 
nonzero travel time spread estimate. In the transition 
region, taken for convenience to be < < A0, one 
may choose a Hamiltonian of the form H = H(I — s0) 
where s = AI/A0, and H = H{I) elsewhere. It follows 
that in the transition region 1(0) = lo + s0 and / = 
constant elsewhere. A simple generalization of the above 
calculation then gives for a complete ray cycle 



-T(r D )« 
(h - lo) 



(2tt - AO) 



g(/oK(/ ) 

M'o)) 2 



A0 
~2~ 



(21) 



and 



p r {r )(r ~ r a ) 
H(Io) 



-(2.-^(1,-/0)^ 



so the sum (recall [T5]l is 



ST^^(h-Io). 



.(22) 



(23) 



It should be emphasized that this expression is first order 
accurate in AI = I\ — Iq, but that no assumption about 
the smallness of A0 has been made. As noted above, 
incomplete cycle end segment pieces give 0((AI) 2 ) cor- 
rections to ST. If A0 has approximately the same value 
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at each upper turn, then one has after n upper turns, 
correct to O(AI), 

A0-A, . A0, 

1=1 

AO 

« — v^(A/) rms . (24) 

Consistent with the numerical simulations shown in Fig. 
12, AO ps 0.8 radians and AI w 4 ms. (We, and in- 
dependently F. Henyey [personal communication], have 
confirmed that these estimates also apply under AET- 
like conditons.) With these numbers and yfn = 8, ap- 
propriate for AET, (|24|l gives a time spread estimate 
of approximately 13 ms, in approximate agreement with 
the numerical simulations shown in Fig. 11. Note that 
(|24|l does not preclude a travel time bias. A cautionary 
remark concerning the use of (|24|l is that Virovlyansky 
has pointed out that, owing to secular growth, the 
0((AI) 2 ) contribution to ST may dominate the O(AI) 
contribution at long range. 

Figure 9 shows that in the AET environment simulated 
near-axial ray time spreads are greater than simulated 
steep ray time spreads. We do not fully understand this 
behavior. A possible explanation for this behavior is that 
time spreads increase as rays become increasingly flat ow- 
ing to the breakdown of the apex approximation. But one 
would expect that this trend should be offset, in part or 
whole, by the relative smallness of internal-wave-induced 
sound speed perturbations near the sound channel axis. 
In addition, we have seen some evidence that an addi- 
tional factor may be important in the AET environment. 
Namely, we have observed a positive correlation between 
travel time spreads and stability exponents; stability ex- 
ponents in the AET environment are shown in Fig. 5. 
We have chosen not to dwell on flat ray time spreads in 
this paper because there are no AET measurements of 
these spreads to which simulations can be compared. It 
is clear, however, that the issues just raised need to be 
better understood. 



IV. WAVEFIELD INTENSITY STATISTICS 

In this section we consider the statistical distribution 
of the intensities of both the early and late AET arrivals. 
Recall that experimentally the early arrival intensities 
have been shown 0, E3 to approximately fit a lognor- 
mal probability density function (PDF) and the late ar- 
rival intensities have been shown |25j to fit an exponen- 
tial distribution. The late arrival exponential distribu- 
tion is not surprising as this distribution is characteristic 
of saturated statistics. The early arrival near-lognormal 
distribution is surprising, however, inasmuch as theory 
|3 EH1 predicts saturated statistics, i.e., an exponential 
intensity PDF. It has been argued [3^| that the theory 
can be modified in such a way as to move the early ar- 
rival prediction from saturated to unsaturated statistics. 



The latter regime is characterized by a lognormal inten- 
sity PDF. This fix is conceptually problematic inasmuch 
as, in this theory, the unsaturated regime is character- 
ized by the absence of micromultipaths, which seriously 
conflicts with the numerical simulations presented in the 
previous section where the number of micromultipaths 
is very large. In this section we provide self-consistent 
explanations for both the late arrival exponential distri- 
bution and the early arrival lognormal distribution. The 
challenge is to reconcile the early arrival near-lognormal 
intensity PDF with the presence of a large number of 
micromultipaths. Some of the arguments presented are 
heuristic, and some build on the results of numerical sim- 
ulations. A complete theoretical understanding of inten- 
sity statistics has proven difficult . 

Our approach to describing wavefield intensity statis- 
tics builds on the semiclassical construction described by 
Eq. (J2J. At a fixed location it is seen that the wave- 
field amplitude distribution is determined by the distri- 
bution of ray amplitudes and their relative phases. Note 
that both travel times and Maslov indices influence the 
phases of ray arrivals. Complexities associated with tran- 
sient wavefields and caustic corrections will be discussed 
below. 

An important observation is that in the AET envi- 
ronment, including internal- wave-induced sound speed 
perturbations, simulated geometric amplitudes of both 
steep and flat ray arrivals approximately fit lognormal 
PDFs. This is shown in Fig. 14. Previously it has 
been shown [2(j that ray intensities in a very different 
chaotic system also fit a lognormal PDF; in that system 
single scale isotropic fluctuations are superimposed on 
a homogeneous background. (Note that all powers of a 
lognormally distributed variable are also lognormally dis- 
tributed, so ray amplitudes have this property if and only 
if ray intensities have this property.) The apparent gen- 
erality of the near-lognormal ray intensity PDF suggests 
that it applies generally to ray systems that are far from 
integrable; the arguments presented in Ref. [2(| suggest 
that this should be the case. The intensity distributions 
of the early and late AET arrivals, corresponding to steep 
and flat rays, respectively, will be discussed separately. 

We consider first the early arrivals. We saw in the pre- 
vious section that the micromultipaths that make up one 
of these arrivals have the same ray identifier and have 
a very small spread in travel time. The dominant mi- 
cromultipaths were seen (see Fig. 11) to have a time 
spread of approximately 1 ms which is a small fraction 
of the approximately 13 ms period of the 75 Hz carrier 
wave. Also, our simulations show that the Maslov in- 
dices of the dominant micromultipaths differ by no more 
than one unit. These conditions dictate that interference 
among the dominant micromultipaths is predominantly 
constructive. Because travel time differences are so small, 
the pulse shape should have negligible influence on the 
distribution of peak intensities. In a model of the early 
AET arrivals consisting of a superposition of interfer- 
ing micromultipaths, the micromultipath properties that 
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play a critical role in controlling peak wavefield intensi- 
ties are thus: 1) their amplitudes have a near-lognormal 
distribution; 2) the dominant micromultipaths have neg- 
ligible travel time differences; and 3) the dominant mi- 
cromultipaths have Maslov indices that differ by no more 
than one unit, ft should be noted that very different 
behavior would have been observed if: the source band- 
width were significantly more narrow as this would have 
caused micromultipaths with different ray identifiers to 
interfere with one another; or the source center frequency 
were significantly higher as phase differences between in- 
terfering micromultipaths would then have been signifi- 
cant. 

An additional subtlety must be introduced now: the 
PDF of the intensities of the constituent micromulti- 
paths that make up a single arrival is not identical to 
the PDF described in Ref. |2Cj and shown in the upper 
and middle panels of Fig. 14. The latter PDF describes 
the distribution of intensities of randomly (with uni- 
form probability) selected rays leaving the source within 
some small angular band and whose range is fixed. This 
PDF is biased in the sense that it overcounts the micro- 
multipaths with large intensities and undercounts those 
with small intensities. Unbiased micromultipath inten- 
sity PDFs can be constructed from the biased PDFs 
shown. To do so, consider a manifold (a smooth curve in 
phase space corresponding to a fan of initial rays) which 
begins near (zq,pq) and arrives in the neighborhood of 
(z r ,p r ) at range r. Sampling in fixed steps of the differ- 
ential Spo (as was done to produce the upper and mid- 
dle panels of Fig. 14; this is equivalent to uniform ran- 
dom sampling) leads to a highly nonuniform density of 
points on the final manifold since (zo,po + Spa) prop- 
agates to (z r + q2iSpo,p r + qnSpo). The greater q 2 i, 
the lower the density of points locally at final range. 
A uniform sampling in final position is acheived instead 
by considering the initial conditions that would lead to 
(z r + Sz r ,p r + qn/q2iSz r ). Its density of points on the ini- 
tial manifold can be deduced from its initial condition, 
(zq,Po + Sz r /q2i). It is necessary to sample (721 times 
more densely on the initial manifold in order to acheive 
uniform sampling in Sz r at range r. To account for this 
effect, we need to know the PDF for (?2i with uniform 
initial sampling. Roughly speaking, the PDF of the ab- 
solute values of the individual matrix elements of q have 
the same form as for \Tr(Q)\, apart from a shift of the 
centroid that is lower order in range than the leading 
term. From the results of Ref. 20], the (biased in the 
sense described above) probability that 521 falls in the 
interval between x and x + dx is 
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1 
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2r{D - v L ) 



, x > . (25) 



value of g2i) and vl is the true Lyapunov exponent. The 
new (unbiased in the sense described above) PDF, p', for 
uniform sampling in Sz r is related to the previous one by 



p'\ q2 i\( x ) 



;P\g 2 i\( x ) 



2nr(i> — vl) x 



exp 



(In (x) — vr) z 
2r{v-v L ) 



, x > (26) 



Here v is a finite-range estimate of the Lyapunov ex- 
ponent based on an average (over an ensemble of rays) 



where the factor x accounts for the extra counting weight 
of Q2i, and (x) just preserves normalization and is calcu- 
lated using Eq. (|25(l . This calculation shows that the 
unbiased micromultipath ray intensity PDF also has a 
lognormal distribution; the only change relative to the 
biased PDF is an increase in the mean from vi,r to Dr. 
Because lognormality is maintained, this correction rep- 
resents only a trivial change to the problem. 

Approximate lognormality of the constrained (eigen- 
ray) PDF of ray intensity is shown in the lower panel 
of Fig. 14. This PDF was constructed using the same 
eigenrays that were used to produce Figs. 8 and 11. The 
corresponding unconstrained ray intensity PDF is shown 
in the middle panel of Fig. 14. (Two constraints - fixed 
receiver depth and fixed ray identifier - are built into 
the lower panel PDF. A very similar constrained PDF 
results if only the receiver depth constraint is applied, 
provided ray launch angles are limited to the 'steep' ray 
band used to construct the middle panel.) It should be 
noted that the constrained (eigenray) PDF shown in the 
lower panel of Fig. 14 was constructed from numerically 
found eigenrays; because weak eigenrays are difficult to 
find numerically they are undercounted and the PDF is 
biased. In a practical sense this bias is of little conse- 
quence because the weak eigenrays that are difficult to 
find contribute negligibly to the wavefield. 

We return now to the problem of simulating the early 
AET arrivals. It is tempting to think that because the 
constituent micromultipaths have a near-lognormal dis- 
tribution, the sum of their contributions should also be 
near-lognormally distributed. Unfortunately, this is gen- 
erally not the case. Consider, for example, the special 
case in which phase, including Maslov index, differences 
are negligible. Then all micromultipaths interfere con- 
structively and peak wavefield amplitudes can be mod- 
elled as the sum of many lognormally distributed vari- 
ables. Because all moments of the lognormal distribu- 
tion are finite, the central limit theorem applies. Under 
these conditions, if sufficiently many contributions are 
summed, the distribution of the sums - the wavefield 
amplitudes - would be a gaussian. 

To simulate the statistics of the early AET arrivals (re- 
call Fig. 14 and the accompanying discussion) we have 
used several variations of a simple model. An arrival 
was modelled as a sum of n m interfering micromultipaths 
whose: 1) amplitudes are lognormally distributed; and 
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2) phases, uTi — jiiTv/2 mod 2w, have a clearly identifi- 
able peak. Micromultipath contributions were coherently 
added. The peak intensity of the sum - whose travel time 
is not known a priori - was then recorded. Using an en- 
semble of 10 4 peak intensities, a peak intensity PDF was 
then constructed. Peak intensity PDFs constructed in 
this fashion were found to be very close to lognormal; a 
typical example is shown in Fig. 15. In this example 
/ = lo/2-k = 75 Hz, the TVs were identical, \ii € j,j + 1 
with equal probability (note that choice of the integer j 
is unimportant) and n m — 5. Other combinations of dis- 
tributions for Ti (either a Gaussian or the limiting case 
of a delta distribution), [ii (taken either from j,j + 1 or 
j — 1 , j, j + 1 with equal probability) , and the parameter 
n m (between 2 and 100) were tested. These simulations 
showed that provided the phase constraint noted above 
was satisfied, a near log-normal peak intensity PDF re- 
sulted. 

Two points regarding this simple model are notewor- 
thy. First, this model does not constitute a theory of 
wavefield peak intensity statistics, but it does serve to 
demonstrate that our ray-based simulations of the early 
AET arrivals are consistent with the measured distribu- 
tion of peak intensities. Second, simulations (not shown), 
performed with E = 1.0 Eqm yield sets of dominant mi- 
cromultipaths that violate assumption 2); phases are uni- 
formly distributed and summing micromultipath contri- 
butions yields a distribution of peak intensities that is 
not close to lognormal. Thus, our simulations suggest 
that a near-lognormal distribution of early arrival peak 
intensity requires a relatively weak internal wave field. 

A complication not accounted for in the preceeding dis- 
cussion is the presence of caustics. At caustics geometric 
amplitudes JHJ diverge and diffractive corrections must 
be applied. At short range (on the order of the first focal 
distance - a few tens of km in deep ocean environments) 
we expect that intensity fluctuations will be dominated 
by diffractive effects. The entire wavefield should be or- 
ganized by certain high-order caustics which leads to a 
PDF of wavefield intensity with long tails [39l |40j, |4l| . In 
spite of the importance of diffractive effects at short range 
- and probably also at very long range - we believe that, 
in the transitional regime described above, intensity fluc- 
tuations are not dominated by diffractive effects. This 
somewhat counterintuitive behavior can be understood 
by noting that in the vicinity of caustics the importance 
of diffractive corrections to JSJ decreases as the curva- 
ture of the caustic increases. Under chaotic conditions 
the curvature of caustics increases, on average, with in- 
creasing range, so the fraction of the total number of 
multipaths that require caustic corrections decreases, on 
average, with increasing range. This is true even as the 
number of caustics grows exponentially, on average, in 
range. This argument leads to the somewhat paradoxi- 
cal conclusion that, prior to saturation at least, we expect 
that the importance of caustic corrections decreases with 
increasing range. 

We turn our attention now to the late unresolved AET 



arrivals, corresponding to the near-axial rays. Here, time 
spreads are sufficiently large that micromultipaths with 
different ray identifiers are not temporally resolved, i.e., 
are not separated in time by more than (A/) -1 . At each 
(z,T) in the tail of the arrival pattern the wavefield can 
be modelled as a superposition of micromultipath con- 
tributions with random phases. The quadrature compo- 
nents of the wavefield have the form of sums of terms of 
the form 

xi = a,iCOs(4>i) , and 

Ui = a i sin(</) i ) (27) 

where (f>i is a random variable uniformly distributed on 
[0, 2tt). The distribution of ai is close to lognormal, but a 
correction must be applied to account for pulse shape as 
many of the interfering micromultipaths partially over- 
lap in time. This correction is unimportant inasmuch as 
the central limit theorem guarantees that, provided the 
distributions of Xi and yi have finite moments, the dis- 
tributions of sums of Xi and yt converge to zero mean 
gaussians. Thus, wavefield intensity is expected to have 
an exponential distribution, consistent with the observa- 
tions. The comments made earlier about caustics apply 
here as well. 

The question of what causes the transition from the 
structured early portion of the AET arrival pattern to 
the unstructured finale region deserves further discussion. 
Recall that the early resolved arrivals have small time 
spreads and peak intensities that are near-lognormally 
distributed, while the finale region is characterized by 
unresolved arrivals and near-exponentially distributed in- 
tensities. In both regions the time spreads and intensity 
statistics are consistent with each other inasmuch as in 
our simulations a near-lognormal intensity distribution 
is obtained only when there is a preferred phase, while 
the exponential distribution is generated when phases are 
random, i.e., when phases are uniformly distributed on 
[0,2tt). 

With these comments in mind, it is evident that the 
most important factor in causing the transition to the 
finale region is the increase in internal-wave-scattering- 
induced time spreads as rays become less steep; as 
time spreads increase, neighboring timefront branches 
blend together and the phases of interfering micromulti- 
paths get randomized. The trend toward increasing time 
spreads as rays become less steep is evident in Fig. 9. 
The surprising result is that the scattering-induced time 
spreads of the steep arrivals are so small; we have seen 
that this can be explained by making use of the apex ap- 
proximation. As noted at the end of the previous section, 
we do not fully understand the cause of the trend toward 
larger time spreads as rays flatten. In the finale region 
internal-wave-scattering-induced time spreads exceed the 
time difference between neighboring timefront branches 
that would have been observed in the absence of internal 
waves. Fig. 1 shows that these time gaps decrease as 
rays become increasingly flat. Indeed, this figure shows 
that even in the absence of internal waves there would 
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not have been any resolvable timefront branches in the 
last half second or so of the arrival pattern. Some loss 
of temporal resolution in the measurements is of course 
due to the finite source bandwidth, but without internal 
waves (or some other type of ocean fluctuations) wave- 
field phases would not be random; instead a stable inter- 
ference pattern would be observed. 

Finally, we note that an interesting feature of our ray 
simulations is that the near-axial rays have much higher 
stability exponents (typically about (100km) -1 ) than the 
steeper rays (typically about (300 km) -1 ); see Fig. 5. 
Also, note that Figs. 6 and 7 strongly suggest that the 
near-axial rays in the AET environment are much more 
chaotic than the steeper rays. The relative lack of sta- 
bility of the near-axial rays in the AET environment is, 
we believe, caused by the background sound speed struc- 
ture. This topic will be discussed in detail elsewhere. 
We have chosen not to focus on this topic here becuse 
the measured intensity statistics in the AET finale re- 
gion are not very sensitive to the near-axial ray intensity 
PDF; as noted above, the argument leading to the expec- 
tation that wavcficld intensity in the finale region should 
have an exponential distribution holds for a very large 
class of ray intensity distributions. 

Finally, we note that an interesting feature of our 
ray simulations is that the near-axial rays have much 
higher stability exponents (typically about (100 km) -1 ) 
than the steeper rays (typically about (300 km) -1 ); see 
Fig. 5. Also, note that Figs. 6 and 7 strongly sug- 
gest that the near-axial rays in the AET environment 
are much more chaotic than the steeper rays. The rel- 
ative lack of stability of the near-axial rays in the AET 
environment is, we believe, caused by the background 
sound speed structure. This topic is discussed in de- 
tail in Ref. [iij]: there it is shown that ray stability 
is largely controlled by a property of the background 
(which is assumed to be range-independent) sound speed 
profile, a — (I / uj)<Lj / dl . Large values of \a\ are asso- 
ciated with ray instability, a vs axial ray angle in the 
range-averaged AET environment, and in five 650 km 
block range-averaged sections of the AET environment, 
are shown in Fig. 16. (The choice of averaging over 650 
km blocks in range is, of course, arbitrary. Range aver- 
aging should be done locally however, because the local 
sound speed structure may be very different than that 
which results after averaging over the entire propagation 
path.) The relatively large near-axial ray values of |a| 
seen in this figure are consistent with the strongly chaotic 
nature of these rays seen in Figs. 5, 6 and 7. Unfortu- 
nately, the measured wavefield intensity statistics in the 
AET finale region are not very sensitive to the near-axial 
ray intensity PDF; as noted above, the argument leading 
to the expectation that wavefield intensity in the finale 
region should have an exponential distribution holds for 
a very large class of ray intensity distributions. Thus, we 
are not aware of any way that the AET measurements can 
be used to test our finding that the near-axial rays have 
larger stability exponents, on average, than the steeper 



rays. 

In summary, we believe that the observed near- 
lognormal PDF of wavefield intensity for the early re- 
solved AET arrivals is transitional between fluctuations 
dominated by caustics at short range and saturation at 
long range, where phase differences will be larger. Sur- 
prisingly, for the steep ray AET arrivals phase differences 
among the steep ray AET arrivals are very small so that 
saturation has not been reached. The late-arriving AET 
energy, on the other hand, is characterized by interfer- 
ing micromultipaths with random phases, leading to an 
exponential PDF of wavefield intensity. Here, the under- 
lying near-lognormal PDF of ray intensity is obscured. 
We believe that diffractive effects do not control the in- 
tensity fluctuations of either the early or late AET ar- 
rivals. Although our theoretical understanding of many 
of the issues raised in this section is clearly incomplete, it 
is worth emphasizing that our ray-based numerical sim- 
ulations of intensity statistics, in which ray trajectories 
are predominantly chaotic, are in good qualitative agree- 
ment with the AET measurements. We do not believe 
that this agreement is accidental. 



V. DISCUSSION AND SUMMARY 

In this paper we have seen that in the AET environ- 
ment, including internal-wave-induced sound speed per- 
turbations, ray trajectories are predominantly chaotic. In 
spite of extensive ray chaos, many features of ray-based 
wavefield predictions were shown to be both stable and 
in good agreement with the AET measurements. Pre- 
dicted and measured spreads of acoustic energy in time, 
depth and angle were shown to be in good agreement. It 
was shown that associated with the chaotic motion of ray 
trajectories is extensive micromultipathing, and that the 
micromultipathing process is highly nonlocal; the many 
micromultipaths that add at the receiver to produce what 
appears to be a single arrival may sample the ocean very 
differently. It was shown, surprisingly, that on the early 
arrival branches the highly nonlocal micromultipathing 
process causes only very small time spreads and does not 
lead to a mixing of ray identifiers. A quantitative ex- 
planation for the cause of the very small time spreads of 
the early ray arrivals was presented. Partially heuristic 
explanations for the near-lognormal and exponential dis- 
tributions of wavefield intensities for the early and late 
arrivals, respectively, were provided. 

The fact that one is able to make apparently robust 
and accurate predictions of many wavefield features us- 
ing ray methods under conditions in which ray trajecto- 
ries are predominantly chaotic may surprise some read- 
ers. Chaotic motion is, after all, intimately linked with 
unpredictability. This apparent paradox is reconciled by 
noting that while individual chaotic ray trajectories are 
unpredictable beyond some short predictability horizon, 
distributions of chaotic ray trajectories may be quite sta- 
ble and have robust properties |45j . Consider, for exam- 
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pie, the motion of the ray whose initial conditions are 
z(0) = -652 m, (p(0) = 8.5°, in the AET environment in- 
cluding a known internal-wave-induced sound speed per- 
turbation. At r — 3252 km, the depth z, angle <p, travel 
time T, and even the number of turns M made by this 
ray are likely to be, for all practical purposes, unpredicat- 
blc. In contrast, the distribution in (z, ip, T, M) of a large 
ensemble of rays with z(0) = -652 m, 8° < <p(0) < 9° is 
very stable in the sense that essentially the same distri- 
bution is seen for any large ensemble of randomly chosen 
(with uniform probability) rays inside this initial angu- 
lar band. It is the latter property that allows us to make 
meaningful predictions using an ensemble of chaotic rays 
in a particular realization of the internal wave field. In 
addition, ray distributions with similar statistical proper- 
ties are observed using different realizations of the inter- 
nal wave field, suggesting that these statistical properties 
of rays are robust. 

Our exploitation of results that relate to ray dynam- 
ics, including ray chaos, leads to a blurring of the dis- 
tinction that is traditionally made between deterministic 
and stochastic wave propagation problems. The ray dy- 
namics approach emphasizes the distinction between in- 
tegrable and nonintegrable ray systems, corresponding to 
range-independent and range-dependent environments, 
respectively. In generic range-dependent environments 
at least some ray trajectories will exhibit chaotic motion. 
The oceanographic origin - mesoscale variability, internal 
waves or something else - of the range-dependent struc- 
ture is not critical. An important conceptual insight that 
can come only from exploitation of results relating to ray 
chaos is that generically phase space is partitioned into 
chaotic and nonchaotic regions. This mixture contributes 
to a combination of limited determinism and constrained 
stochasticity. 

For those who wish to exploit elements of determinism 
in the propagation physics for the purpose of perform- 
ing deterministic tomographic inverses, the results that 
we have presented have important implications. We have 
seen, for example, that, in spite of extensive ray chaos, 
many ray-based wavefield descriptors are stable and pre- 
dictable, and should be invertible. A less encouraging but 
important observation is that although the travel times 
of the steep early arrivals are stable and can be inverted, 
the nonlocality of the micromultipathing process limits 
one's ability to invert for range-dependent ocean struc- 
ture. 

For those who wish to understand and predict wave- 
field statistics, our results also have important implica- 
tions. These comments are based on our analysis of the 
AET measurements, but are expected to apply to a large 
class of long-range propagation problems. First, we have 
seen that stable and unstable ray trajectories coexist and 
that this influences wavefield intensity statistics. Sec- 
ond, we have seen that micromultipathing is extensive 
and highly nonlocal. The strongly nonlocal nature of the 
micromultipathing process is significant because: a) this 
process cannot be modelled using a perturbation analy- 



sis which assumes the existance of an isolated background 
ray path; and b) the effects of this process will not, in 
general, be eliminated as a result of local smoothing by 
finite frequency effects. Third, we have seen that the 
appearance of stochastic effects is not entirely due to in- 
ternal waves. We have seen, for example, that, even in 
the absence of internal waves, near- axial rays in the AET 
environment are chaotic. Recall that we noted in the 
introduction that part of the motivation for the present 
study was to understand the reasons underlying the find- 
ing of Colosi et al. [24| that the theory described in Ref. 
|31| fails to correctly predict the AET wavefield statistics. 
This failure is not surprising in view of the observation 
that this theory is based on assumptions that either ex- 
plicitly violate, or lead to the violation of, all three of 
the aforementioned properties of the scattered wavefield. 
Note also in this regard that the combination of extensive 
micromultipathing, small time spreads, and lognormally 
distributed intensities that characterizes the early AET 
arrivals is not consistent with any of the three propaga- 
tion regimes identified in this theory. 

Although we have argued that a ray-based wavefield 
description which relies heavily on results relating to ray 
chaos can account for all of the important features of 
the AET measurements, it should be emphasized that 
our analysis has some shortcomings. Recall that our 
simulations using E = 0.5 Eqm resulted in steep ar- 
rival time spreads and peak intensity statistics in good 
agreement with measurements, but that simulations us- 
ing E = 1.0 Eg m did not, although direct measurements 
favor E = 1.0 Eqm- We have not provided a rigorous 
argument to establish the connection between a near- 
lognormal ray intensity PDF and a near-lognormal wave- 
field intensity PDF. Indeed, a complete theory of wave- 
field statistics which accounts for complexities associated 
with ray chaos is lacking. This is closely linked to the sub- 
ject of wave chaos [2^ , which is currently not well under- 
stood. Also, we have only briefly discussed the spreads 
of energy in depth and angle, and more work needs to be 
done on quantifying time spreads. 

Finally, we wish to remark that the importance of ray 
methods is not diminished by recent advances, both the- 
oretical and computational, in the development of full 
wave models, such as those based on parabolic equa- 
tions. The latter are indispensible computational tools 
in many applications. In contrast, the principal virtue 
of ray methods is that they provide insight into the un- 
derlying wave propagation physics that is difficult, if not 
impossible, to obtain by any other means. The results 
presented in this paper illustrate this statement. 
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FIG. 1: Measured and simulated AET wavefields in the time-depth plane. Upper panel: a typical measurement, shown on a 
gray scale plot with a dynamic range of 30 dB, of wavefield intensity. Middle panel: ray simulation with internal waves. Lower 
panel: ray simulation without internal waves. Wavefield intensity in the simulations is approximately proportional to the local 
density of dots, each corresponding to a ray. 
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FIG. 2: Same as Fig. 1 except that plane wave beamformed wavefields are shown. The beamformed ray simulations assume 
a dense receiving array whose upper and lower bounds coincide with the bounds of the AET receiving array. The reference 
depth was taken to be the depth of the uppermost hydrophone, z = —0.9 km. 



17 



tfic [m a A ] 




FIG. 3: Measured and simulated rms sound speed fluctuations in the AET environment. The measurements shown were made 
near the source (labelled east) and near the receiving array (labelled west). The simulated fluctuations are due to internal 
waves using both E = 0.5 Eqm and E = 1.0 Eqm, as described in the text. 
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FIG. 4: Simulated Lagrangian manifolds in the AET environment without (left panel) and with (right panel) internal waves. 
The ray density used in the plot on the right and for small \p\ in the plot on the left is too sparse to resolve what should, in 
each plot, be a smooth unbroken curve that does not intersect itself. 
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FIG. 8: Three segments of the range-depth plane along the AET transmission path showing a subset of a set of micromultipaths. 
The three rays that have the strongest amplitudes are shown using dashed lines; these rays are not resolved on this plot. 
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FIG. 9: Ray depth vs. travel time in the AET environment with internal waves, but only for rays with identifiers +137 (left 
panel) and +151 (right panel). The points plotted are a subset of those plotted in the lower panel of Fig. 1. 
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FIG. 10: Ray identifier vs. launch angle for a subset of the rays that are plotted in Fig. 6. 
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FIG. 11: Upper panel: normalized ray intensity vs. travel time for a subset of the eigenrays (micromultipaths) with identifier 
+137 that connect the AET source and a receiver at depth 1005 m at the AET range. Middle panel: the corresponding Maslov 
indices vs. travel time. Lower panel: envelope of the waveform synthesized by coherently adding the contributions from all the 
micromultipaths shown (solid curve), and the normalized envelope of the AET source waveform (as seen in the far field; dashed 
line) . The three rays whose intensity is largest are shown using solid circles in the upper and middle panels; the remaining rays 
are shown using open circles. 




FIG. 13: Schematic diagram showing a twice apex-scattered ray (heavy line) and the path that the same ray would have 
followed in the absence of apex scattering events (light line) . 
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FIG. 14: Cumulative probability density of ray intensity in the AET environment. Upper panel: flat rays (\ip\ < 5°) sampled 
uniformly in launch angle. Middle panel: steep rays (6° < \ip\ < 11°) sampled uniformly in launch angle. Lower panel: eigenrays 
with +137 identifier connecting the AET source and a receiver at depth 1005 m at the AET range. In all three panels the 
dashed lines correspond to lognormal distributions whose mean and variance approximately match those of the simulations. 
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FIG. 15: Cumulative probability density of peak wavefield intensity simulated using a simple model as described in the text. 
The same simulated cumulative density is compared to exponential (upper panel dashed line) and lognormal (lower panel 
dashed line) cumulative density functions. This figure should be compared to Fig. 14 in Ref. |24| . 
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FIG. 16: Stability parameter a vs. axial ray angle in the range- averaged AET environment (dashed curve) and in five 650 km 
block range-averages of the AET environment (solid curves). 



